Import des packages

rm(list=ls())
library('dplyr')
library('tidyr')
library('VennDiagram')
library('grid')
library('futile.logger')
library("ggplot2")
library("gplots")
library('plotly')
library('tools')
library('cowplot')
library('UpSetR')

Import des VCF

Import des 3 vcf obtenus après lancement du pipeline (vcf_12x, vcf_24x, vcf40x) et import du vcf de référence (vcf_ref).

vcf_12x = read.table("/home/elodie/Documents/Analysis/gold_12x_on_data_elodie_filtervarianttranches_2D_decomposed_normalized_uniq.vcf")
vcf_24x = read.table("/home/elodie/Documents/Analysis/gold_24x_on_data_elodie_filtervarianttranches_2D_decomposed_normalized_uniq.vcf")
vcf_40x = read.table("/home/elodie/Documents/Analysis/gold_40x_on_data_elodie_filtervarianttranches_2D_decomposed_normalized_uniq.vcf")
vcf_ref = read.table("/home/elodie/Documents/Elodie/Datas_ismael_ref/gold_on_data_elodie_ismael_decomposed_normalize_uniq.vcf")

Exploration des VCF et des variants

Nombre et type de variants par VCF

plot1 = ggplot(nombre_variants, aes(fill=type, y=vcf, x=nombre)) +
  geom_bar(position = "stack", stat = "identity") + 
  scale_fill_manual(values = c("#bebada", "#ffffb3", "#8dd3c7", "#170f0f")) +
  labs(fill = "Type de variants", x = "Nombre de variants", y = "VCF") +
  theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5)) + 
  labs(title = "Nombre et type de variant par VCF") 
ggplotly(plot1)

La distribution des types de variants semble équivalente pour chaque VCF. Cependant dans les 3 VCF générés par le pipeline on trouve des variants “Undetermined” mais pas dans le VCF de référence.

Nombre de variants par gène et par VCF

plot2 = ggplot(nombre_variants_par_gene, aes(fill=gene, y=vcf, x=nombre)) +
  geom_bar(position="stack", stat="identity") +
  scale_fill_manual(values = c("#fdc086", "#91bfdb", "#b2df8a")) +
  labs(fill = "Nom des gènes", x = "Nombre de variants", y = "VCF") + 
  theme(panel.grid = element_blank(),plot.title = element_text(hjust = 0.5)) + 
  labs(title = "Nombre de variants par gène et par VCF") 
ggplotly(plot2)

Le nombre de variants est identique dans chaque VCF pour les gènes UTS2 et LPL mais le nombre varie pour le plus grand gène EPHA5.

Nombre et type de variants par gène et par VCF

nombre_variants_par_gene$vcf_nom_genes = paste(nombre_variants_par_gene$vcf, nombre_variants_par_gene$gene, sep="_")

plot3 = ggplot(nombre_variants_par_gene, aes(fill=type, y=vcf_nom_genes, x=nombre)) +
  geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
  scale_fill_manual(values = c("#bebada", "#ffffb3", "#8dd3c7", "#170f0f")) +
  labs(title = "Nombre de variants par gène, par VCF et par type de variant") +
  labs(fill = "Type de variants", x = "Nombre de variants", y = "VCF") + 
  theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))
  
ggplotly(plot3)

La répartition des types de variants par gène et par VCF est concordante. Aucune insertion dans UTS2 dans les 4 VCF.

Exploration des SNV

Nombre et type de SNV par VCF

plot4 = ggplot(nombre_type_snv[!is.na(nombre_type_snv$type), ], aes(fill = type, y = vcf, x = nombre))+
  geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
  scale_fill_manual(values = c("#fc9272", "#9ebcda")) +
  labs(title = "Type de SNV par VCF", x = "Type de substitutions", y = "Nombre de substitutions") +
  labs(fill = "Type de SNV") +
  theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))
plot4

En génétique humaine, il y a un environ 2 fois plus de transitions que de transversions. Sur le plot, on voit environ 2/3 de transitions et 1/3 de transversions ce qui est cohérent. Valeurs cohérentes dans tous les VCF.

Nombre et type de transitions/transversions par VCF et par gène

plot5 = ggplot(nombre_type_snv[!is.na(nombre_type_snv$type_snv), ], aes(fill = type_snv, y = vcf_gene, x = nombre)) +
  geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
  labs(title = "Répartition des types de SNV par VCF et par gène", x = "Nombre", y = "VCF et gène") +
  scale_fill_manual(values = c("#a6cee3", "#1f78b4","#b2df8a", "#33a02c", "#fb9a99", "#e31a1c", "#fdbf6f", "#ff7f00", "#cab2d6", "#6a3d9a", "#ffff99", "#b15928")) +
  labs(fill = "Type de SNV", x = "Nombre", y = "VCF") +
  theme(plot.title = element_text(hjust = 0.5, size=14, face="bold"), plot.subtitle = element_text(hjust = 0.5)) + 
  theme(panel.grid = element_blank())
ggplotly(plot5)

Une répartition des types de SNV cohérente entre VCF pour chaque gène.

Exploration des délétions/insertions

plot6 = ggplot(vcf_all %>% filter (type == "insertion" | type == "deletion"), aes(x = taille_delins, y = vcf)) + 
  geom_boxplot(aes(fill=type)) + 
  theme_classic() + 
  labs(x = "Taille des delins", y = "VCF", title = "Boxplots de la taille des insertions et délétions par VCF") + 
  theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) + 
  labs(fill = "Type de delins") +
  scale_fill_manual(values = c("deletion" = "#bebada", "insertion" ="#ffffb3"))
plot6

La distribution des tailles des délétions est équivalente entre les 4 VCF. Les VCF 12x, 24x et 40x contiennent plus de valeurs aberrantes pour les insertions et délétions que le VCF de référence. La taille maximale des insertions du VCF de référence est plus petite que pour les 3 VCF obtenus avec le pipeline. L’interquartile est plus petit également pour la taille des insertions du VCF de référence ce qui signifie que les valeurs des insertions sont plus regroupées autour de la médiane. Il y a moins de disparité des valeurs pour le VCF de référence.

plot7 = ggplot(vcf_all %>% filter (type == "insertion"), aes(x = taille_delins, y = vcf)) +
  geom_boxplot(aes(fill=gene)) +
  theme_classic() +
  labs(x = "Taille des insertions", y = "VCF", title = "Boxplots de la taille des insertions par VCF et par gène") + 
  theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) + 
  labs(fill = "Noms des gènes") +
  scale_fill_manual(values = c("EPHA5" = "#fdc086", "LPL" ="#91bfdb", "UTS2" = "#b2df8a"))

plot8 = ggplot(vcf_all %>% filter (type == "deletion"), aes(x = taille_delins, y = vcf)) +
  geom_boxplot(aes(fill=gene)) +
  theme_classic() +
  labs(x = "Taille des délétions", y = "VCF", title = "Boxplots de la taille des délétions par VCF et par gène") + 
  theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) + 
  labs(fill = "Noms des gènes") +
  scale_fill_manual(values = c("EPHA5" = "#fdc086", "LPL" ="#91bfdb", "UTS2" = "#b2df8a"))

plot_grid(plot7, plot8)

Exploration de la profondeur de lecture des variants

plot9 = ggplot(vcf_12x_24x_40x, aes(x = as.numeric(vcf_12x_24x_40x$DP), y = vcf)) + geom_boxplot(aes(fill=gene)) + theme_classic() + 
  theme_classic() +
  labs(x = "Profondeur de lecture des variants", y = "VCF", title = "Boxplots de la profondeur de lecture des variants par gène") + 
  theme(plot.title = element_text(hjust = 0.5, size=12, face="bold"), plot.subtitle = element_text(hjust = 0.5)) + 
  labs(fill = "Noms des gènes") +
  scale_fill_manual(values = c("UTS2" = "#b2df8a", "LPL" ="#91bfdb", "EPHA5" = "#fdc086"))
plot9
## Warning: Use of `vcf_12x_24x_40x$DP` is discouraged.
## ℹ Use `DP` instead.

plot10 = ggplot(dp_12x, aes(x = DP, y = nombre)) +
  geom_histogram(stat = "identity", fill = "#a6cee3") +
  geom_vline(xintercept = mean(dp_12x$DP), color = "#fd8d3c", linetype = "dashed", size = 1) +
  labs(x = "Profondeur de lecture", y = "Nombre de variants", title = "VCF 12x") + 
  theme(axis.text.x = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5, size = 14, face = "bold")) + 
  theme(panel.grid = element_blank())
## Warning in geom_histogram(stat = "identity", fill = "#a6cee3"): Ignoring
## unknown parameters: `binwidth`, `bins`, and `pad`
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
plot11 = ggplot(dp_24x, aes(x = DP, y = nombre)) +
  geom_histogram(stat = "identity", fill = "#1d91c0") +
  geom_vline(xintercept = mean(dp_24x$DP), color = "#fd8d3c", linetype = "dashed", size = 1) +
  labs(x = "Profondeur de lecture", y = "Nombre de variants", title = "VCF 24x") + 
  theme(axis.text.x = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5, size = 14, face = "bold")) + 
  theme(panel.grid = element_blank())
## Warning in geom_histogram(stat = "identity", fill = "#1d91c0"): Ignoring
## unknown parameters: `binwidth`, `bins`, and `pad`
plot12 = ggplot(dp_40x, aes(x = DP, y = nombre)) +
  geom_histogram(stat = "identity", fill = "#253494") +
  geom_vline(xintercept = mean(dp_40x$DP), color = "#fd8d3c", linetype = "dashed", size = 1) +
  labs(x = "Profondeur de lecture", y = "Nombre de variants", title = "VCF 40x") + 
  theme(axis.text.x = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5) ) + 
  theme(panel.grid = element_blank())
## Warning in geom_histogram(stat = "identity", fill = "#253494"): Ignoring
## unknown parameters: `binwidth`, `bins`, and `pad`
plot_grid(plot10, plot11, plot12, ncol = 3, nrow = 1)

Exploration de la VAF des VCF 12x, 24x et 40x

On ne fera pas d’exploration de la VAF ni de la profondeur DP pour le VCF de référence car il n’y a pas d’explication aux valeurs de DP, AD, ADALL dans le VCF. Les valeurs ne sont pas cohérentes avec ce qui est observé dans le bam de référence associé au VCF dans IGV.

VAF minimale détectée dans chaque VCF

cat("VAF minimale dans VCF 12x :", min(vaf_12x$vaf), "%", "\n")
## VAF minimale dans VCF 12x : 16.7 %
cat("VAF minimale dans VCF 24x :", min(vaf_24x$vaf), "%", "\n")
## VAF minimale dans VCF 24x : 14.3 %
cat("VAF minimale dans VCF 40x :", min(vaf_40x$vaf), "%", "\n")
## VAF minimale dans VCF 40x : 16.3 %
plot13 = ggplot(vcf_vaf, aes(x = nombre_variants, y = cat_vaf , fill = vcf)) +
  geom_bar(position = position_fill(reverse=FALSE), stat="identity") +
  labs(title = "Répartition du nombre de variants par catégorie de VAF", x = "Nombre de variants", y = "Catégorie de VAF") +
  scale_fill_manual(values = c("#a6cee3", "#1d91c0", "#253494")) +
  labs(fill = "VCF") +
  theme(plot.title = element_text(hjust = 0.5, size=14, face="bold"), plot.subtitle = element_text(hjust = 0.5)) + 
  theme(panel.grid = element_blank()) +
  theme_classic()
ggplotly(plot13)

Exploration des variants en commun dans les 4 VCF

variants_list = list(vcf_12x = unique(vcf_12x$variants), vcf_24x = unique(vcf_24x$variants), vcf_40x = unique(vcf_40x$variants), vcf_ref = unique(vcf_ref$variants))

plot14 = upset(fromList(variants_list), order.by = "freq", text.scale = 2,  point.size = 5, sets.bar.color = "#999999", main.bar.color = "#fc8d62")
plot14

plot15 = display_venn(variants_list, category.names = c("vcf_12x" , "vcf_24x" , "vcf_40x", "vcf_ref"), 
             fill = c("#a6cee3", "#1d91c0", "#253494", "#E69F00"), cat.cex = 1.2,
        cat.fontface = "bold",
        cat.default.pos = "outer",
        cat.dist = c(0.1, 0.1, 0.1, 0.1))

plot16 = ggplot(cat_vcf_all_nombre, aes(fill=type, x=nombre, y=categorie)) +
  geom_bar(position = "stack", stat="identity") +
  scale_fill_manual(values = c("#bebada", "#ffffb3", "#8dd3c7", "#170f0f")) +
   labs(title = "Nombre de variants par catégories et par type de variants") +
  labs(fill = "Type de variants", x = "Nombre de variants", y = "Catégories de VCF") + 
  theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))
plot17 = ggplot(cat_vcf_all_genes, aes(fill=gene, x=nombre, y=categorie)) +
  geom_bar(position = "stack", stat="identity") +
  scale_fill_manual(values = c("UTS2" = "#b2df8a", "LPL" ="#91bfdb", "EPHA5" = "#fdc086")) +
  labs(title = "Nombre de variants par catégories et par gènes") +
  labs(fill = "Type de variants", x = "Nombre de variants", y = "Noms des gènes") + 
  theme(panel.grid = element_blank(), plot.title = element_text(hjust = 0.5))

plot_grid(plot16, plot17, ncol = 1, nrow = 2)

Exploration des variants en commun dans les 3 VCF générés par le pipeline

variants_list_12x_24x_40x = list(vcf_12x = unique(vcf_12x$variants), vcf_24x = unique(vcf_24x$variants), vcf_40x = unique(vcf_40x$variants))

plot18 = display_venn(variants_list_12x_24x_40x, category.names = c("vcf_12x" , "vcf_40x", "vcf_24x"),
  fill = c("#a6cee3", "#1d91c0", "#253494"), cat.cex = 1.2,
        cat.fontface = "bold",
        cat.default.pos = "outer",
        cat.dist = c(0.1, 0.1, 0.1))

Exploration des métriques de qualité: datas générés par Hap.py d’Illumina

recall (sensibilité), précision (spécificité) et score F1. Ce sont des métriques qui permettent de voir si un modèle est performant. Le score F1 combien recall et précision Plus le recall est élevé, plus le modèle repère des positifs. Plus la précision est élevée, moins le modèle se trompe sur les positifs (moins il y a de faux positifs). Plus le score F1 est élevé, plus le modèle est performant.

Courbe ROC

happy_12x_summury = read.table("/home/elodie/Documents/Analysis/Happy_12x/Happy_12x/ref-12x.summary.csv", header=TRUE, sep=',')
happy_24x_summury = read.table("/home/elodie/Documents/Analysis/Happy_24x/Happy_24x/ref-24x.summary.csv", header=TRUE, sep=',')
happy_40x_summury = read.table("/home/elodie/Documents/Analysis/Happy_40x/Happy_40x/ref-40x.summary.csv", header=TRUE, sep=',')
plot_roc_snp = plot_data(roc_snp, is.PR=TRUE)
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
plot_roc_indel = plot_data(roc_indel, is.PR=TRUE)
plot_grid(plot_roc_snp, plot_roc_indel, nrow=2, ncol=1)